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Diffraction of Sound by a Half-Plane in a Uniform Flow 

R. K. Amiet 


SUMMARY 

The classical problem of the diffraction of sound by a half-plane is extended to the case where the half-plane is 
immersed in a uniform flow. Both the leading edge and trailing edge cases axe considered. The sound source is assumed 
to be a point source. Computer programs are presented for the evaluation of the necessary integrals. The results are 
expressed relative to both the sound that would be present in the absence of the half-plane and the sound at the point on 
the edge where the ray passes the edge. A minor misrepresentation of the solution in the literature has been clarified for 
the case t r > 0, which corresponds to the illuminated region. Several sample calculations are presented. 
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INTRODUCTION 


This document is Volume II of a 5 volume report on aerodynamics and noise of advanced turboprops. Volumes 
I* IE and IV relate to aerodynamic and acoustic disturbances caused by the blades and Volume V relates to shielding of 
propeller noise in the cabin by the fuselage boundary layer. This volume presents theory and a computer code developed 
for evaluation of the shielding benefit that might be expected by an airplane wing in a wing-mounted propeller 
installation. Several computed directivity patterns are presented to demonstrate the theory. Application to actual airplane 
geometry is given in Volume HI. 

The diffraction of a sound or electromagnetic wave by a half-plane is a classic problem in the literature. The 
first rigorous soludon for the case of a plane wave incident on a half-plane was given by Sommerfeid (ref. 1). Other 
authors (refs. 2-6) have subsequently treated this problem. This fundamental soludon has been extended by many authors 
to other cases such as the diffraction of a point source field by a half-plane, the diffraction of a plane wave by a wedge, 
the case where the half-plane is considered "soft** (the boundary condition of no-flow through the surface is replaced by 
some other boundary condition), etc. 

The primary interest of these solutions has been for electromagnetic applications, with acoustics being a 
secondary interest. For the electromagnetic problem there is no such thing as a mean flow, so the case of diffraction in a 
flow has received little attention. Thus, it has only been recently with the advent of the concept of using the wing of an 
aircraft for noise shielding that the case of diffraction by a surface in a flow has been given attention. 

The present analysis is based on the case of diffraction with no flow. By combining a Galilean and a Lorentz 
transform, the wave equation with a mean flow can be reduced to the ordinary wave equation. The boundary conditions 
must also be transformed, but because of the simple geometry of the problem the boundary conditions remain unchanged 
in the transformed plane. A similar transformation was used by Candel (ref. 7). 

Allowance is also made in the analysis for the case of a swept wing. The same combination of Galilean and a 
Lorentz transforms mentioned above lead to a problem with no flow but a different sweep. Thus, with proper 
interpretation and reverse transformation of the zero flow case, one can obtain the solution for the effect of sweep. This 
is included in the computer program. 

The solution procedures for the cases of leading and trailing edges are basically the same. For the leading edge 
one works with the velocity potential since the boundary condition is on normal velocity. For the trailing edge the 
boundary condition in the wake is that there be no pressure jump. For this case the same solution applies, but the 
fundamental variable is the pressure. 

Two normalizations of the solution are given by the computer program. The first is just to normalize the 
solution by the pressure that would exist at the observer in the absence of the plate. The second choice is to normalize 
the solution by the pressure at the point on the edge at which the ray crosses in passing to the observer. This edge 
pressure calculated in the absence of the plate is extrapolated to be an equal distance from the source as the observer. The 
reason for this second normalization is to better represent those cases where the point source is very directional. The 
pressure at the observer in the absence of the plate bears no relation to the diffracted field at the observer in the presence 
of the plate, whereas the pressure at the edge crossing point is on the ray that moves from the source to the observer. 

FORTRAN computer programs are presented with detailed documentation. The output from these programs 
compares favorably with the results of other investigators. 
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ANALYSIS 


A. Zero Mean Flow 


1. Exact solution 

The solution for the sound field of a source near a semi-infinite plane without flow has been given previously 
(refs. 1 - 6) so that here all that need be done is to give a summary of the results. The geometry of the problem is shown 
in figure 1. Mathematically the problem can be stated as follows. 

A source of strength Q, exp(i©t) at (x 0 ,y 0 ,zj induces a velocity potential <D 0 that satisfies the wave equation 

(v 2 + * 2 )d> 0 = Q 0 S(x-xJ 8(y-yj5(z-zj 0) 

where k = ©/c 0 . The plate is assumed to lie in the y = 0 plane in the region x > 0. The boundary condition of no flow 
through the plate is 

v y = d&fjdy = 0 ony = 0 f x>0 (2) 


The solution to this problem can be written (ref. 4) 


<p =Q 0 iL H™(kR cosh fi)dn+ I H?\kS cosh fi) dp 

87C IU L 


(3) 


where 


R = V/f + r 2 - 2r 0 r, cos(6 0 - 6,) + (z 0 - zf S ■ + rj - 2v, c<«(0 o + 6} + (z 0 - 7f 

hr H, - SM' [i (4±4j 


The source coordinates are x, = (x,,y,,zj (or (r.,0,,z,) in polar coordinates) and the observer coordinates are x o = (x o .y o ,zJ 
or (rA.zJ. 

As a check on the solution, consider placing the observation point near the source; i.e., x 0 — » x,. Then R — > 0 
and m -♦ oo. The second integral in equation (3) remains bounded while the first goes to infinity as R 1 . Because of the 
fae.mr R in the argument of the Hankel function, it is not immediately obvious that the upper limit p* of the integral can 
be set equal to infinity in the R -» 0 limit. It can be shown, however, that the major contribution to the integral comes 
from the region p * 0; using the relation 

-i //,® (a coshu) dn = !-«■“ (4) 

2 a 


( see Bom and Wolf (ref. 6, p.586)) equation (3) becomes 

4> 0 “» - GcA 4 **) + G(l) as R -> 0 (5) 

Then, near the source, v= a<P 0 /9R = Q^itR 2 ). Considering a small sphere around the source and multiplying the 
velocity by the area 4nR*, the total volume flow out of the sphere is found to be Q,. Thus, equation (3) represents the 
diffraction field of a monopole source of strength near a half plane. This is consistent with equation (1), the right 
hand side of which can be shown to represent a monopole. 
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2. Integral approximations 


The integrals in equation (3) can be approximated for large k. The two integrals are similar in form so that only 
the approximation for the first will be discussed. This integral will be analyzed in the two regimes (^>0 and p^ < 0. 
For p, < 0 both limits on the integral are positive. For p* > 0, the integral will first be divided into two integrals; one 
with limits i°° and one with both limits negative; the integral with infinite limits can be evaluated in closed form, and 
the integral with both limits negative is treated in the same way as for the case when both limits are positive. 

Thus, first consider the case p* < 0 for which both limits are positive. It follows from the definitions of p, and 
R in equation (3) that 

R cosh n R = */(.r 0 + rf + (z o - z f > R , (6) 

Thus, over the entire range of integration (for which the following inequality holds: 

kr cosh pi>kR x ( 7 ) 

For kRj»l the Hankel function can be replaced by its asymptotic form for large argument giving 


I = I H®\kR cosh pi) dpi-- 



V cosh pi 


where the parameter x, defined as 



e 




Y(t 2 +1)(t 2 + 2) 


dt 


T = VJ sink (p /2) 


( 8 ) 


(9) 


was introduced in place of p. The definition of p, from equation (3) gives for the corresponding value of X, 


*r = 2 


V 


r , r . 


R{R + /?.) 


cos 



V /?,//? - 1 

- fRJR - 1 


e a - e,<x 
a 0 -e^7[ 


( 10 ) 


Equation (8) can be approximated further. Since kR is large, variation with p of the exponent kR cosh p is also 
large, and the method of stationary phase can be used to evaluate the integral. Since the derivative of the phao. i s zero 
only at the point p = 0 and since this point is not included within the range of integration, there are no stationary phase 
points. Following Senior (ref. 4), x is set equal to its lower limit in the non exponential part of the integral. The 
justification for this can be found in the book by Erdelyi (ref. 8), p.51, under a discussion of the method of stationary 
phase. For integrals with rapidly varying phase the major contribution comes from the end points and any stationary 
phase points. Because the portion of the integrand under the radical in equation (8) varies slowly compared to the phase of 
the exponential, it is justifiable to set the argument under the radical equal to the constant value taken at the end point of 
the integration. Substitution from equation (10) in equation (8) gives the approximate value valid for large kR, 


where 


/- 


2 e -»(«♦«*) 
kR 


R F \-z R i2kR/7C ) 

Y#,(*,+*) 


T ,<0 


(ID 


F (*) = 

F is related to the Fresnel integrals C(x) and S(x) by 


e-^dt 


(12a) 


F'<?) = i-d.-E\x) 


(12b) 
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where 


E'(x) h C W - i 5 (x) = 



(12c) 


Equation (1 1) is valid only for t R < 0. The case t R > 0 is treated simply by dividing the integration range into 
two parts as discussed previously. The integral with limits ± <*> is found in closed form using equation (4). The integral 
with both limits negative is treated in the same manner as for positive limits, and is easily found to give the same result. 
The final result for t r > 0 is 


/ <• e"*** -2_ 

kR 


i + e" 


iRyiR.+R) 


F' {xJlkRjn) 


T „>0 


(13) 


For x R = 0 Eqs. (1 1) and (13) become identical. The first term in equation (13) is the field of a source in free space, while 
the second gives the contribution of the edge. For a physical interpretation of the difference between Eqs. (11) and (13), 
note that from equation (10) for 0 O - 0, < n that equation (13) applies when 0 O - 0, < n; i.e. t when the observer is not in 
the shadow of the plate, while equation (1 1) is for an observer in the shadow. For the second term in equation (3), an 
equation similar to equation ( 13) applies for an observer that can see the reflection of the source in the plate (0 o + 0. < 
jc), and an equation similar to equation (1 1) applies for an observer who cannot (0 O + 0 § > re). 

The author has not found in the literature the two different forms which were given here for the two regions x R > 
0 and x R < 0. Thus, Senior (ref. 4, eq. (33)) and Bom and Wolf (ref. 6, eq. (45) of section 1 1.7) give the same result as 
equation (11) here, but do not make a distinction for the case x R > 0. In Eqs. (1 1) and (13) above, the function F* always 
has a positive argument. The results given by Senior and by Bom and Wolf have a negative argument for F when x R > 
0; this can be rewritten with a positive argument by dividing the range in the integration of F* in equation (12a) into two 
ranges, one with limits ±°° and the other with both limits negative (which is the same as both limits positive). The 
function F* will then be the same as in equation (13) above, but the integrated part will have a factor R/fR 1 (R 1 + R)/2] w 
that does not appear in equation (13). This may not have been noticed earlier since for positive x R the discrepancy 
between equation (11) and equation (13) is most apparent in the iluminated region where the diffraction effects are 
miminftl and so not of interest to calculate. Near the diffraction zone x R = 0 and Eq. (10) shows that R, - R so that the 
factor R/fR,^, + R)/2] lfl is nearly equal to 1. 

Both the exact and the approximate solutions are antisymmetrical about H* = x R = 0. Thus, only the region x R 
< 0 need be considered in detail. If the results are plotted against the variable X R (k*R) 1/2 then the argument of F in 
equation (1 1) is given by the abscissa. If the integral is multiplied by the factors appearing before F* in equation (11), 
the asymptotic form of the integral will be independent of k*R when plotted; the deviation of the exact solution from this 
approximate curve will then show more clearly than if I in Eqs. (8) and (1 1) is plotted directly. This plot given in figure 
2 clearly shows that for values of k-R > 1 the approximate result is very close to the exact integral. 


B. With Mean Flow 


1. Unswept wing 

The addition of a mean flow to the problem will now be considered. Without flow, and for a general time 
dependence the problem can be formulated in terms of the velocity potential as the following boundary value problem: 

(' - -V ^r) f* = Q <' > 5 <* - **) 5 o - yJ s ( 2 * O ( 14 > 

l 4 dl 2 l 


The boundary conditions are 


30 /dy = 0 on y = 0, x > 0 


(15) 
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For the problem with flow the edge of the plate will be assumed to be a leading edge, as opposed to a trailing 
edge which will be discussed later. The coordinate system is fixed to the plate with the source being stationary and the 
fluid moving toward increasing x; see figure 1. The equation for the velocity potential <> becomes 


-M-+ 

df 


M 


dx 


4>~ QU) 8(x- x t ) 8(y -y t ) 8(z - zj 


(16) 


which can easily be verified by performing a Galilean transformation fixing the observer to the fluid in which case the 
left hand side reduces to equation (14). The boundary conditions remain as before in equation (15). 


The equation can now be transformed using the following combination of the Lorentz and Galilean 
transformations: 

x^X, y-+Y/l l 

where p 2 m 1 - M 2 . Equation (16) becomes 


z -* Z //}, t -* T - M X /(co/J 2 ) 


& 


V„ 2 - -L 


fa 37*1 


i = Q(T- k 'Mx} 8 (x - x,) 8(Y /fj-yj 6 (Z //) - z t ) 


(17) 


(18) 


where V 0 is the gradient operator in the transformed coordinates X,Y,Z and 

<t> 0 (XJZ,T ) = 0Cr f y f z t / ) k* = k Ip 2 (19) 

For any constant c the following property of generalized functions holds true: 

S(cx) = c’ 1 S(x) (20) 

Applying this equation and assuming a sinusoidal time dependence 

e(O = Q 0 ^ f 0o(O=<*>o^ (21) 

gives 

(?o + k* 2 ) 0 O = e, 8 (X - X} 8 (Y - YJ 8 (Z - Z,) (22a) 

where 

Cl = ^2oe‘ a * A£,, (226) 

Equation (22a) for M * 0 is identical in form with equation (1) for the M = 0 case. The amplitude of the source 
strength Q, which is a constant in equation (1) becomes Q, in equation (22a); Q, is also a constant in that it does not 
depend on the independent variables X,Y,Z. Also, the boundary conditions in equation (15) are unchanged by the 
transformation. Denote the solution for the case M = 0 by < t , 0 (x,,y t ,z i ;x <i ,y <> ,z o ,k). The solution of equation (16) for the 
forcing function given by equation (21) is thus found from equation (3) by replacing Q„ by Q„ exp(-ik’Mx^, k by k* and 
x,y,z by X,YJZ respectively. Finally the result is multiplied by exp(io)T) = expflcot + ik*Mx). After inverse 
transformation from X.Y.Z to x,y,z the result is 

<Kx J .y,rZ,^ a .y 0 .z 0 ,k,M,t) = **> <J> 0 (.x^pz^py^k) e im (23) 

A similar transformation was used by Candel (ref. 7). His transformation and final result can be seen to be 
essentially the same as the above. For example, his k, = k/fi and x, = x/P; thus, his expfi^Mx,) is equivalent to 
exp(ik MxJ above. However, because of the factor k that appears outside the square brackets in equation (3), and 
the definition of k* used here differs from that of Candel, a comparison with Candel's equation (28) might appear to differ 
from that obtained above. The difference results from the fact that Candel is considering the case of an incident plane 
wave which is assumed to be the same with and without flow while the present case is concerned with the case of a point 
source with an output which is affected by the flow. 
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2. Swept wing 


The potential field of a swept wing in a flow satisfies the same equation, equation (16), as for the unswept case. 
The only difference is in the boundary condition. Equation (15) becomes 

30 /3y = 0 on y = 0, x - z tan y> 0 (24) 

where y represents the angle between the leading edge and the z axis; see figure 3. The same transformation given by 
equation (17) then reduces the problem to the no flow case. The sweep angle y ' in the transformed zero-flow plane is 

tan y'= /T 1 tan y (25) 

The solution 0 for the potential field of the source near a swept wing with flow can then be written 

<&x s% y a j^ 0 d 0 a 0 \kM*yt ) = &'o(x J &,fo J * 0 'Py 0 foo> k 'y) * icot (26) 

where <P 0 ’ represents the solution to the no-flow problem in a coordinate system with nonzero y 9 ; note the prime on <D 0 ' 
to distinguish it from O which was defined in a coordinate system with Z along the leading edge. <t> does not have y* as a 
parameter. The fact that y ' is nonzero means that in order to relate to the previous solution given by equation (3), a 
coordinate rotation through the angle y* must be made. The transformed coordinates, denoted by primes, are related to the 
non-primed coordinates in the zero-flow plane by 


X’ = X cos y f - Z sin y' Z’ = X sin y'+ Z cos y’ 


(27) 


The equation relating the velocity potential for flow to that with no flow and with the z axis along the leading edge is 
finally 


<KXsWx oM ;kju 9 yj) = 0 o (x;,y>X,y>> 


(28) 


3. Leading and trailing edge cases 

Equation (29) gives a solution to the wave equation which satisfies the boundary condition of no flow through 
the plate. For the case of a leading edge (0 < M < 1) the solution represents the velocity potential, not the pressure. The 
pressure field for the leading edge case can be found from the momentum relation 

p = - p D<f>fDt (29) 

where D/Dt = 3/3t + U 3/3x represents the substantial derivative. The velocity potential for the leading edge case must 
behave as x w , in contrast to the pressure and velocity fields which have a x‘ w singularity there. Thus, there is no 
discontinuity in potential at the leading edge, but there is a discontinuity in pressure. The behavior of the velocity 
potential at the leading edge is typical for this type of problem. The same behavior is found, for example, for a 
two-dimensional plate moving along its normal in an incompressible fluid; this can readily be solved using complex 
variable mapping techniques. 

For the trailing edge case, the behavior is different. Here the imposition of the Kutta condition requires that the 
pressure field go to zero at the trailing edge. Also, there can be no pressure discontinuity downstream of the trailing 
edge, although there will, in general, be a discontinuity in the velocity potential. Looking again at the boundary value 
problem for the leading edge given by Eqs.(15) and (16), the solution is seen to apply to the trailing edge if M — > - M 
and if <j> is taken to represent the pressure field. Thus, with -1 < M < 0, equation (28) represents the pressure field of a 
point pressure monopole in the presence of a trailing edge. With 0 < M < 1, equation (28) represents the velocity 
potential of a point velocity-potential-monopole near the leading edge of a semi-infinite zero-thickness flat plate. 

Comparing the solutions for the leading and trailing edges, it will be noted that they are solutions for different 
forcing functions. The leading edge solution is for a monopole in the velocity potential. The trailing edge solution is 
for a monopole in the pressure field. The monopole in the velocity potential can be interpreted as fluid injected into the 


7 



scream through a porous body. The monopole in the pressure field can be interpreted as the same type of source injecting 
fluid through a porous body but with the addition of a dipole aligned with the flow; another interpretation is that it 
represents a body in the fluid which pulsates in thickness but which adds no fluid to the flow, in contrast to the first case 
of fluid added through a porous body. 

Because of the difference in source type for leading and trailing edges, it is necessary to normalize the solutions 
with respect to the incident pressure field, in some manner, before a comparison between the leading and trailing edge 
solutions can be made. 


4. Normalization of the solution 

In order to properly gauge the importance of diffraction using the preceding solution, it is necessary to first 
normalize the solution in some manner. As it stands, equation (3) predicts a diffracted amplitude which is directly 
proportional to the source strength Qq. The most obvious method of normalizing the solutions for leading and trailing 
edge diffracted pressure is to divide by the pressure that would be present at the observer location in the absence of the 
half plane. With flow, the solution to the monopole forcing function is 

= _ Op £()) 

where 

= V(*o - xf + J 1 [(y 0 - yf + (z 0 - z/ 


For the trailing edge case 4> 0 be taken to be the pressure while for the leading edge case the substantial derivative must 
be taken. In taking the derivative, only the far-field component will be retained. It is felt that the additional complexity 
introduced by retaining the near-field terms in the normalizing factors would obscure the meaning of the results. Both the 
near and far-field terms are present in the actual solution, however, and it is a simple matter to account for the near-field 
terms in the normalizing factor if desired. 

For the trailing edge case the normalizing factor is 


For the leading edge case 


with 




observer ' 


Qq 

4X0 r o , 


N U=Po 


£ >± 

Dt 


•Pi o 


observer 


(1 -McosO.) 

4no ot p 2 c 0 


COS 0J = xJO" 


(31) 

(32) 


These factors should be used to normalize the appropriate diffraction relations given previously. The computer 
program, to be described later, makes this normalization using the exact diffraction solution. (Actually, if kR, is large, 
the asymptotic solution is used, but this loses little in accuracy.) 

Before settling on this normalization for the analysis, let us first consider further the properties of the diffraction 
solution. For present purposes, to simplify matters, let us consider just the asymptotic solution for kRj large and for M 
= 0. Then, from Eq.(3) and Eqs.(ll) - (13) 

*=G 0 (*/8ff) [/(*) + /($)] (33) 


I(R)-e M 2L 

kR 

.«J*i-** 


1 -e 


Re-**'-* 

2i2*kr/Jt,cos[{6 0 -e,m 


_L 


kf^A i2nkr s cas[{e 0 -8,y 2] 


**>0 

**<0 


(34) 
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with the restriction that 0 O - 0 § not be near 7t; i.e., the observer must not be near the edge of the diffraction field. The 
diffraction component of this result is 


l/R) 



4V?i + i*4 


Combining this with I d (S) gives 

^ Qo l cos ( 6J2 ) cos (6J2) 

( 2 n) m Ikr/J ^ cos + cos Q s 


(35) 

(36) 


For M = 0, because of the omnidirectionality of the source, it makes no difference whether the normalization is with 
respect to the pressure at the observer position or at the edge (extrapolated to an equal distance). Also, the solutions for 
the leading and trailing edges become identical. If equation (36) is normalized using equation (32) with M = 0, the result 


cos {BJ2) cos (6/2) 


cos Q 0 + cos 6 S 


This represents the diffraction solution for either the leading or trailing edge case (assuming M = 0) normalized by either 
the pressure at the observer with no plate or the pressure at the edge (extrapolated to equal distance from the source) with 
no plate. It represents both the normalized potential solution and the normalized pressure since as noted from equation 
(29) the only difference between the two is the constant factor -icop which is present in both the solution and the 
normalizing factor, and so drops out. Because it is only the diffraction component of the solution, the direct and reflected 
rays must be added if the observer is not assumed in the shadow zone of the diffraction. 

This result is very interesting. Consider the directional behavior of this normalized field. For an observer in the 
plane z Q = z # the angular dependence occurs only in the explicit 0 O factor in equation (37). Also, recall that there is no 
directionality in the normalizing term. The directionality in equation (37) is the same as that for the case of a plane wave 
incident on an edge. (See, e.g., equation (13) of the paper by Candel (ref. 7), remembering to account for the fact that 
Ca nde l defines 0 and 0 to be the supplements of 0, and 0 § defined here.) If the observer goes to the far field, r 0 » r >t the 
amplitude factor (multiplying the directivity) becomes [2r B /(rck)] l/2 /r 0 . and this can readily be shown to be the same for all 
forcing functions such as dipole, quadrupole, etc.; i.e., take the derivative of the monopole solution and normalize as 
above with the incident pressure field. 

The fact that the same directivity pattern is obtained from two completely different source types (monopole and 
plane wave) is illustrative of a very powerful result: for the limits of high frequency and far field 

kR x » 0 and r 0 » r t (38) 


the diffraction field for an arbitrary source directivity can be found from that for a monopole source. The two diffraction 
patterns are equal if normalized by the wave amplitude incident on the edge; that is, the diffraction pattern is dependent 
only on the ray striking the edge, and not on the general source directivity. The relevant edge point is the one which 
minimizes the dis ta nce from source to edge to observer. Further discussion of this principle can be found in the paper by 
Keller (ref. 9). It is also described in the report by Boeing (ref. 10). 

This is basically a geometric acoustics type of solution in which the incident sound moves along a ray tube to a 
point on the edge. On reaching the edge, which is assumed thin with respect to the wavelength, the sound diffracts 
around the edge. Now, the diffraction from any point on the edge is determined by the properties of the incident wave 
within approximately one wavelength of the point; under the assumption kr B » 1, the radius of curvature of the incident 
wave is much greater than a wavelength and the wavefront can be treated as plane. Thus, the diffracted component of the 
solution for the limits given by equation (38) is independent of the detailed source directivity, only depending on the 
amplitude directed at the edge. Once reaching the edge, the ray is diffracted into a conical surface with the z axis being the 
axis of the cone. For any point on this conical surface, the minimum distance from this point to the source, with the 
edge as an intermediate point, is along a path from this point to the apex of the cone, and from the apex of the cone to 
the source. If the observer is in the shadow of the plate, then the only sound reaching the observer is the diffracted sound; 
if the observer is not in the shadow zone of the plate, then the direct sound (and the reflected sound if the observer is in 
the shadow zone of the image source) must be added to the diffracted sound. In this case, of course, the source directivity 
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at these angles, in addition to that of the ray striking the edge, are important. 

For nonzero M, the same results apply. The ray paths in the transformed plane correspond with those in the 
original plane. Since the solutions for a plane wave and a point source agree for M = 0, they will agree for M * 0 since 
the solution for nonzero M is found by transforming both the plane wave solution and the monopole source solution for 
M = 0 using the same transformation. 

Since it is the ray impinging on the edge which determines the diffraction solution, it makes more sense to 
normalize the solution in terms of the source pressure (in the absence of the plate) at the edge of the plate rather than that 
at the observer. This will allow the development of a solution for a general source directivity (under the assumptions in 
equation (38)) to be simply obtained from that for the monopole solution. Since we wish to calculate the sound 
reduction produced by introducing the plate, it might be thought that this would be given by the solution normalized by 
the pressure at the observer without the plate. This is true only for the monopole source since this is the source type for 
which the solution was derived. If a more general solution is desired, the source intensity at the edge must be used for 
normalization. From this solution normalized by the sound intensity at the edge, the attenuation of the sound from a 
particular source type produced by introducing a plate can be calculated if the ratio of the directivities at the edge and at 
the observer are known for this general source type. Also note that if the pressure at the observer is used for the 
normalization, the normalization pressure would depend on observer location. 

In order to make this normalization, first the point at which the ray strikes the edge must be determined. 
Denoting this point by (0,0,0* the total path length travelled by the ray from source to observer with this edge point as 
an intermediate point is 

L = Jxf+ y} + (z, - C ) 2 + V x] + y} + (z a - £ ) 2 (39) 

This distance is minimized by setting the derivative 9L/d£ = 0 giving 

C SB ( r o z , + z SW r o + r } (4°) 

For a specific value of £ there is a linear relation between r o and z o . Thus, the locus of observer points related to a specific 
value of £ form a cone; that is, a ray will travel from the source to a point £ on the edge, and from the edge it will 
diffract into a conical pattern. 

This relation is for the zero-flow case without sweep so that for the case with flow and sweep the equation 
would need to be transformed by Eqs.(18) and (26). Denote this transformed position on the edge by (x e ,0,z e ). The ray 
from the source toward this point is 

A = -i (*,-*,) -jy s -k (z,-z,) (41) 

where i , j , k are unit vectors along x, y, z. The vector h x = r o h/lhl is normalized to a length r o to give a vector from the 
source which grazes the edge and ends in the far field. This normalization is to avoid introducing an r' 1 effect on the 
solution due to the difference in distances from the source to the edge and to the observer. With this extrapolation of the 
vector h the normalized far-field pressure for a monopole source with M = 0 will directly give the attenuation produced by 
introducing the plate. For more general sources or for M * 0, the ratio of the directivity in the directions of source and 
edge must be given to get this result. The pressure at this far-field point, to be used for normalization of the result, is 
from equation (32) for a leading edge 

N, = 

4 nc B a, 

where 

= a «l 0 . frtA*) = V (X, - X,) 2 + 0 2 \yf + (z, - zf] (43) 

For the trailing edge, from equation (31) 



N «... = -r 2 - ± V (x s - xf + yj + (z, - zf 
4no. R 


( 44 ) 
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In summary, the pressure at an arbitrary field point due to a monopole at an arbitrary source point is obtained 
from the solution for zero flow using the transformations described above. This is then normalized by the monopole 
pressure at a distance r o from the source along the line touching the edge (assumed to be in the far field). This must then 
be multiplied by the ratio of pressure along the line touching the edge to that at the observer, but for the particular source 
being considered. 


COMPARISON WITH OTHER SOLUTIONS 


Since several of the results presented here have been presented previously by other authors, where possible it is 
worthwhile to compare results with these previous calculations, if for no other purpose than to verify the comp uter 
programs presented here. 


The zero Mach number case shown in figure 4 for the directivity of a point source in the vicinity of a half plane 
can be compared to a similar calculation with the same input parameters in figure 18 of reference 10. The comparison is 
favorable, with complete agreement to the accuracy that the graph can be read. 

The case with flow is slightly more complex. Reference 10 presents figures for leading and trailing edge cases 
with their parameter 0 # = 90*. Here 0 W is used to denote the angle 0 O used in reference 10. From their figure 17 it 
appears that 0 W is the wavefront angle, not the angle to the source position; i.e., 0„ is the included angle between c 0 and 
U. This can be related to the angle 0, drawn to the source position by recalling that the sound ray will propagate at 
velocity c # along its normal and be convected at velocity U along the flow direction x. Let V denote the propagation 
speed of a ray (i.e., the vector addition of the fluid velocity and the acoustic velocity) the following eq uations follow by 
equating the horizontal and vertical components of ray velocity in terms of V and 0 to the values expressed in terms of 
q,, U and 0 W : 

V sin 8 t = Cq sin 8 W V cos 8, - c 0 cos 8 W - U ( 45 ) 

or 

col 8 t = cot 8 W -M /sin 8 W (46) 

The angle 0, is the angle from the edge to the source, and for 0. = 90*. cot 0, = -M. Reference 10 gives results for M = 
if).8; this gives 0, = 141.34* for the leading edge case and 0, = 38.66* for the trailing edge case. 

The results calculated using the present theory for these cases are shown in figures 5a and 5b. The input 
parameters are matched to those of reference 10: k r, = 10, z. = z,, the sweep is zero and the observer is in the far field. 
(Specifically r. = 10 J is assumed in the computer program.) In each of these plots two curves are shown. The only 
difference between the two lines is in the normalization; the dashed curve shows the pressure of the monopole source in 
the presence of the diffracting half-plane normalized by the far-field pressure of the source at the observer with the plane 
removed. The solid curve shows the same diffracted pressure, but normalized by the extrapolated edge pressure as described 
previously. The dashed curves are seen to agree well with the corresponding curves given in figure 18 of reference 10, 
indicating that the results of reference 10 were normalized by the observer pressure with the edge removed. 

As was mentioned previously, the ray striking the edge is the one that determines the diffraction field. Since for 
a general source this ray is independent of the ray toward the observer, it is generally better to normalize the results by 
the ray striking the edge for an observer in the shadow zone. However, for an observer at a location outside the shadow 
zone the ray toward the observer will not be intercepted by the plate; for this case normalization by the ray amplitude in 
the direction of the observer is more meaningful. In fact the curves clearly illustrate that for n-0, <0 <jt + 0 where 
there would be expected to be little effect of the plate (since the observer is outside the shadow zone, but not at a small 
enough angle that he can see the image source), the curve normalized by the observer pressure is in fact near zero, 
whereas this is not so clearly evident from the curve normalized by the edge pressure. Perhaps a combination of these 
no rm a liz at i ons could be made so that when the observer was in the shadow the normalization would be with the edge ray 
while for an observer in the illuminated region the direct ray would be used. 

These curves give the somewhat misleading impression that the edge of the shadow zone gets blown Hark by the 
flow. Recall, however, that the source position for the leading and trailing edge cases is different for the data plotted in 
figures 5. For comparison the results for the case 0 ( 9 90* for both leading and trailing edge are shown in figures 6. 
Here it is evident that the edge of the shadow zone remains at 0 * 270*, unchanged from the M = 0 no flow case. 
Intuitively one might expect the sound to be blown back toward the plate for the leading edge case The explanation of 
this paradox is that the ray propagating toward the edge at the 270* angle must have the normal to its wavefronts canted 
upstream in order for them to propagate at the 270* angle. After passing the edge, the wavefront angle is not changed; 
the ray continues propagation along the same line. 
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BEHAVIOR OF THE SOLUTION 


This section will present sample calculations made using the solution presented previously. Before proceeding, 
however, it is well to determine satisfactory values for the two parameters D and N in the subroutine INTH. D is the 
dividing point between the exact and the approximate solutions and N determines the number of integration points if the 
exact solution is used. 

A value for D can be found from figure 2. The computer program uses the exact or approximate solution based 
on the following relations: 

kR x ^D Exact solution used . 
kR x > D Approximate solution used . 

Using Eq.(10) to replace R, with t*, Eqs.(47) become 

kR (t£ + 1) £ D Exact solution used . 
kR (t^ + 1) > D Approximate solution used . 

The first term in these equations is the abscissa in figure 2. With this fact, the figure indicates that at the points D = 15 
the approximate curve is always very close to the exact curve. This is used as the input to the subroutine INTH. 

The value for N can be determined by doing calculations for increasing N values until no further change takes 
place. With kR = 10, going from a value of N = 5 to N = 10 changed only the fifth digit in the output of the program 
CINTH, which calls INTH for checking the output. The number of integration steps is scaled with kR so that N = 10 
should be a good value for the entire range of kR. The behavior of the solution for variations of the input parameters is 
now considered. 

First consider the variation of the solution with Mach number M. First, compare the leading and trailing edge 
solutions. The only difference in these cases with regard to the computer program is that for the leading edge case M > 0 
is input while for the trailing edge case M < 0 is the appropriate program input. The transformations of x,y,z in 
equation (17) are the same for leading and trailing edge cases. Thus, the basic solutions for the two cases with the same 
value of IM1 have the same magnitude. This is before normalization and before the derivative in equation (29) is taken for 
the leading edge case. If the observer is in the far field, only the phase is important in determining the derivative; the 
same is true for the solution in the absence of the diffracting plane. Also, the variations of the phases for the cases with 
and without the diffracting plane are the same for the same observer location. Thus, when the solutions for pressure are 
normalized by the pressures for an observer at the same point without the diffracting plane, the leading and trailing edge 
cases will give the same results. 

This equality of the solutions for the leading and trailing edges, when normalized by the pressure at the same 
location with no diffracting plane, is evident in figure 7a which is symmetrical about the M = 0 axis. The source is at 0 t 
= 90® and curves are shown for three different observer angles. For simplicity zero sweep and z # = z o are assumed. The 
particular observer angle 0 O = 270° gives the same results for normalization by the edge pressure (extrapolated) and 
observer pressure with no plate since the observer lies on the line from source to edge. Figure 7b shows similar plots 
except for 6 g = 45°. For both these figures there is little dependence on M for the observer on the line joining source to 
edge; i.e., on the edge of the shadow zone. For observer positions further back in the shadow zone there appears to be 
significant dependence on the flow Mach number; generally the effect of M is to increase the shielding for both leading 
and trailing edge cases. 
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THE COMPUTER PROGRAMS 


Included here are two different main programs and several accompanying subroutines. The main program is 
DIFRAC. This program computes the diffraction of a single edge, either leading or trailing edge. This main program is 
used with the subroutines 


Main program 


Main: DIFRAC 

Subroutines: 

SDIFF 

INTO 

INTGC 

BESL 



FC 

CORD 

SRIS 

FRES 


Inputs 

The inputs to the program are defined in terms of the quantities in figures 1 and 3. They are: 

RS = r a = distance of the source from the z axis. 

TS = 0 a = angle of the source from the plane. 

AM = M = Mach number; positive for a leading edge and negative for a trailing edge. 
GAM = y = angle of sweep as shown in figure 3. 

AK = k = co/c 0 = 2 % frequency /sound speed 
RO = r 0 = distance of the observer from the z axis. 

TO = 8 0 =angle of the observer from the plane. 

ZO = z 0 a z coordinate of the observer. 

The inputs can be in any consistent system of units. 


Running the program 

The program begins by prompting the user to enter a name for the file in which the output will be stored, as 
well as being written to the screen. The program then prompts the user to input the variables r f , 0^ M, y, k, r 0 , 0 O , z 0 . 
The program then outputs these variables as a check, and follows with the output of results for either the leading edge 
case if M > 0 or the trailing edge case if M < 0; the pressure, normalized two ways are the first two outputs and the 
following three give the distance of the source from the edge crossing point of the ray. A sample output is 


SOURCE RS = 1.00 THS = 90.0 
2.19322 2.19322 0.0 1.0 

MACH = .000 
.0 

GAM = .0 

AK = 10.0 

RO = 1000.0 

THO = 90.0 

O 

II 

o 

N 

SOURCE RS = 1.00 THS = 90.0 
-5.23460 -5.23460 0.0 1.0 

MACH = .000 
.0 

GAM = .0 

AK = 10.0 

RO = 1000.0 

THO = 270.0 

o 

ll 

O 

N 

SOURCE RS = 1.00 THS = 90.0 
-15.02630 -15.02630 0.0 1.0 

MACH = .000 
.0 

GAM = .0 

AK = 10.0 

RO = 1000.0 

THO = 360.0 

N 

O 

ll 

o 


Compare these with figure 4; the output values 2.19322, -5.23460 and -15.02630 correspond to dB values for the 0 
angles of 90, 270 and 360 degrees. 0 


Program for intermediate verification of theory 

The second main program is CINTH. This is a minor program for checking the subroutine INTH which 
evaluates the integral I in equation (8). This main program is used with the subroutines 

Main: CINTH Subroutines: INTO INTGC BESL 

FC SRIS FRES 
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CONCLUSIONS 


The analytical solution for the diffraction of a plane wave by an infinite span, semi-infinite chord, swept airfoil 
in a mean flow is a straightforward extension of previous results for diffraction without flow or sweep. The results can 
readily be computerized to produce a program that can predict the attenuation of a sound wave produced by wing 
shielding. 
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LIST OF SYMBOLS 


Cp Sound speed 

E f F* Combinations of Fresnel integrals defined by equations (12) 

H t (3> Hankel function of the second kind 

h Vector from the source to point e on the edge; defined in equation (41) 

hj Vector parallel to h with length r o 

I Integral defined by equation (8) 

\ a>/c 0 

k* Transformed value of k; defined in equations (19) 

M Flow Mach number 

P Pressure 

Qff) Source strength as a function of time 

Qo Source strength amplitude; see equations (21) 

Qj Transformed source strength defined by equations (22) 

r 0 Distance of observer from edge 

r s Distance of source from edge 

R Source-observer distance defined in equations (3) 

Rj Distance defined by equation (6) 

S Image-observer distance defined in equations (3) 

v y Fluid velocity normal to plate 

x,y a. Cartesian coordinates defined in figure 1 

X,Y ,Z Transformed coordinates defined by equations (17) 

X',Z' Transformed and rotated coordinates defined by equations (27) 

x,.Zt Point on the edge which gives a minimum for the distance between the source and observer. 
z o z coordinate of observer 

z k z coordinate of source, taken to be zero for simplicity, losing no generality 

P Prandtl-Glauert factor defined in equations (17) 

y Sweep angle defined in figure 3 

Y * Sweep angle in transformed plane defined in equation (25) 

0 O Angle of observer from plane; defined in figure 1 

0 # Angle of source from plane; defined in figure 1 

0 W Used to represent the variable 0 o of Ref. 1 0 

fig, jij Integral limits defined in equations (3) 

p Fluid density 

o M Length parameter defined in equations (30) 

o # Value for o m for an observer on the edge at point e 

t r Parameter defined by equations (9) and (10) 

$ Velocity potential including time dependence 

4> 0 Velocity potential for zero flow including time dependence 

4> d Velocity potential for diffraction component of solution 

<J> 0 Velocity potential for zero flow and with time removed 

co Radian frequency 

Subscripts 

o, s Quantity related to observer or source respectively 
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Appendix A: Main Program for Calculating Diffraction by an Edge 

7/23/84 DIFRAC by R. K. Amiet 


10 

100 


1 
2 

3 

4 

5 

6 

7 

8 

9 

10 
11 
12 

13 300 

14 ] 

15 3 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 


COMPLEX CX, CXI, PHI, PHIP, PHIM, PHIl, PHI1P, PHI1M, 

DPHI , DPHIl , PRLE, PRLE1 , AI 
CHARACTER* 12 FNAME 
AI * CMP LX (0 . , 1 . ) 

WRITE (9,*) f OUTPUT FILE NAME - f 
READ (9,*) FNAME 

OPEN ( 1 , FILE-FNAME , STATUS- • NEW 1 ) 

WRITE (9, 100) 

FORMAT (/' INPUT RS, TS, MACH, GAM, K, RO, TO, ZO, Nl=l END*/) 
READ (9,*) RS, TS, AM, GAM, AK, RO, TO, ZO, N1 
IF (N1 .EQ. 1) GO TO 40 

WRITE (1,300) RS, TS, AM, GAM, AK, RO, TO, ZO 

FORMAT (/ 1 SOURCE 1 , 3X, 1 RS - 1 , F5 . 2, 3X, • THS MACH 

F5 . 3, 3X, ’GAM - ' , F6 . 1, 3X, • AK - 1 , F7 . 1, 3X, 1 RO *\F7.1, 

3X, ’THO - * , F6 . 1, 3X, 1 ZO -\F7.1) 

XS = RS*COS (TS/57 . 2958) 

YS * RS * SIN (TS/57 . 2958) 

ZS - 0. 

B2 - 1.- AM* AM 
B - SQRT (B2) 

AKX = AK/B2 

XO - RO*COS (TO/57.2958) 

YO - RO*SIN (TO/57. 2958) 

CALL CORD (XS, YS, ZS, XO- . 01 , YO, ZO, GAM, B, RSX, TSX, ROX, TOX, ZOX, XE, ZE) 
CALL DIFF1 (RSX, TSX, ROX, TOX, ZOX, AKX, PHIM, PHI1M, 15 . , 10) 

CALL CORD (XS, YS, ZS, XO+ . 01 , YO, ZO, GAM, B, RSX, TSX, ROX, TOX, ZOX, XE, ZE) 
CALL DIFF1 (RSX, TSX, ROX, TOX, ZOX, AKX, PHIP, PHI1P, 15 . , 10) 

CALL CORD (XS, YS, ZS, XO, YO, ZO, GAM, B, RSX, TSX, ROX, TOX, ZOX, XE, ZE) 

CALL DIFF1 (RSX, TSX, ROX, TOX, ZOX, AKX, PHI , PHIl, 15., 10) 

CX - CEXP (AI*AKX*AM* (XO - XS) ) 

CXI - CEXP (AI* AKX* AM*. 01) 

PHI - PHI*CX 
PHIl « PHI 1 *CX 

DPHI - 50 . *CX* (PHIP*CX1 - PHIM/CXl ) 

DPHIl - 50 . *CX* (PHI1P*CX1 - PHI1M/CX1) 

PRLE - - (AK*PHI + AM*DPHI) 

PRLE1 = - (AK*PHI1 + AM*DPHI 1) 

PPTE * 20 . *ALOG10 (CABS (PHI/PHIl) ) 

PPLE - 20 . *ALOG10 (CABS (PRLE/PRLEl) ) 

XOS - XO - XS 

YOS = YO - YS 

ZOS - ZO - ZS 

SGOS • SQRT (XOS*XOS + B2*(YOS*YOS + ZOS*ZOS) ) 

ROS - SQRT ( XOS *XOS + YOS* YOS + ZOS* ZOS) 

XSE » XS - XE 

ZSE - ZS - ZE 

SGSE - SQRT (XSE*XSE + B2* (YS*YS + ZSE*ZSE) ) 

RSE - SQRT (XSE*XSE + YS*YS + ZSE*ZSE) 

CORTE = SGOS*RSE/ (SGSE*ROS) 

CORLE - CORTE* (1 . + AM*XSE/SGSE) / ( 1 . - AM*XOS/SGOS) 

PNORLE - PPLE - 20. *ALOG10 (CORLE) 

PNORTE - PPTE - 20. *ALOG10 (CORTE) 

PP = PPLE 
PNOR = PNORLE 
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55 


IF (AM .LT. 0.) PP - PPTE 


56 


IF (AM .LT. 0.) PNOR - PNORTE 


57 


WRITE (9, 600) PP, PNOR, XSE, YS, 

ZSE 

58 


WRITE (1,600) PP, PNOR, XSE, YS, 

ZSE 

59 

600 

FORMAT ( IX, 2F12 . 5, 3F6. 1) 


60 


GO TO 10 


61 

40 

CLOSE (1) 


62 


END 



COMMENTS: This is the main executive routine for calling the various diffraction subroutines, and for normalizing the 
resulting output by the appropriate pressure. 

3-7 Open file for output. 

8-15 Prompt for inputs. 

RS = Distance of source from z axis. 

TS = Angle of source from plane. 

MACH = Flow Mach no.; M > 0 => leading edge; M < 0 ^ trailing edge. 

GAM = Sweep angle of airfoil; GAM = 0 => z along leading edge. 

AK = co/c 0 

RO = Distance of observer from z axis. 

THO = Angle of source from plane. 

ZO = Distance between z = 0 plane & plane containing observer and normal to z axis. 

16-18 XS, YS, ZS = the x, y, z distances of the source from origin; the origin is located so that the x axis is 
along the flow vector, the y axis is normal to the plane containing the airfoil, and the origin is on the 
leading edge. 

ZS is taken to be zero; this fixes the origin. 

19-20 B = Prandtl-Glauert factor. 

21 AKX = wavenumber in transformed plane. 

22-23 XO &YO = the x and y distances of the observer from the origin. 

24-29 First the transformed coordinates (with an X ending) are found; these are then used in the subroutine 
SDIFF1 for the no-flow solution. The 0.01 deviations from XO are used to find the axial derivative. 

30-31 Phase factor relating the no-flow solution to the flow solution. CXI is the additional factor for solutions 

with 0.01 difference in XO. 

32-33 Relating the flow and the non-flow pressures for the TE case. 

34-35 The x derivatives needed for the calculation of pressure from potential using equation (28). 

36-37 Pressure calculation for the leading edge case. PHI is the pressure for the trailing edge case and PRLE is the 
pressure for the leading edge case. The 1 ending on these two variables denotes the pressure for the case 
with no plate. 

38-39 DB calculations for leading and trailing edge cases after normalizing by the pressure in the absence of the 
plate. 

40-42 Distances between source and observer. 

43 o at observer with source at origin. 

44 Source-observer distance. 

45-46 x and z distances of source from edge-crossing point. 

47 O at edge-crossing point with source as origin. 

48 Distance between edge-crossing point and source. 

49-50 Ratio of pressures at observer and edge-crossing point; line 49 trailing edge case; 

line 50 => leading edge case. 

51-52 DB levels for pressures normalized by that at edge, extrapolated. 

53-56 Print either LE case or TE case, not both. 
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Appendix B: Subroutine for Calculation of Integral in Equation (8) 

8/12/84 INTH by R. K. Amiet 


1 SUBROUTINE INTH (UR, RK, D, N, P, J) 

2 COMPLEX P, AI, E, FTR1 

3 AI = CMPLX (0 . , 1 . ) 

4 RRK = RK*COSH (UR) 

5 IF (RRK .GT. D) GO TO 10 

6 M - N*RRK + 1 

7 CALL INTGC (0 . , UR, RK, M, P) 

8 P = .5* (CEXP (-AI*RK) /RK - AI*P) 

9 J = 0 

10 GO TO 20 

11 10 SQ1 = SQRT ( (RK + RRK) *RRK) 

12 TR = RRK - RK 

13 CALL FRES (TR, E) 

14 FTR1 =» CMPLX ( . 5 , - . 5) - CONGJ(E) 

15 P = FTR1*CEXP (AI* ( .7853982 - RK) ) /SQ1 

16 IF (OR .GT. 0.) P = CEXP (-AI*RK) /RK - P 

17 J = 1 

18 20 RETURN 

19 END 


Comments: This subroutine calculates the integral in the equation following equation (3.2) in the paper by Betty Woods. 

This is the same as the integral I in equation (8) here multiplied by (- i/2). 

I The inputs to the subroutine are UR, RK, D and N. UR = in equation (2), RK = R k with R given by 
equation (3) and k by equation (1). D defines the cutoff point for switching from the exact to the 
asymptotic solution; for k-R t > D, the asymptotic solution is used. The parameter N specifies the number 
of integration points for exact solution; N = 10 appears to be sufficiently large. 

4 RRK = k-R l as can be seen from equation (6). 

5 Check whether to use exact or asymptotic solution. Go to 10 for the asymptotic solution. 

6 M will be proportional to the number of integration points. M increases as N k • R cosh ^ since this is 

proportional to the range of the argument of H t w in equation (8). 

7 The integration subroutine INTGC is called for the integration needed in equation (8). The inputs are the 
integration limits 0 and the parameter Rk and the number M which is proportional to the number of 
integration steps. 

8 Integrated result with limits ± «*> from equation (4) is added to result with limits - ji^.O in equation (8). 

9 J is a program output showing whether the exact or the approximate solution was used. 

J = 0 => exact solution; J = 1 => approximate solution. 

10 Skip over the approximate solution. 

I I SQ1 is the radical [RjfR, + R )] VJ in equation (1 1). 

12-13 Since the subroutine FRES calculates C 2 and S 2 , and since equation (1 1) is defined in terms of C and S, the 
input to FRES must be the argument in equation (11) squared and multiplied by tc/2; see Ref. lip. 300. 
TR used here is t R (R k) w . 

14 FTR1 is F* in equation (11). 

15 The remaining factors in equation (1 1) multiply F\ and the result is multiplied by (- i/2). 

16 For jj* > 0 the final result is composed of an integral from 0 to •• and one from - x R to 0. 

17 Same comment as 9. 
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Appendix C: Subroutine for Calculation of R, S, (x R and m 

8/3/M SDIFF by R. K. Amiet 

1 SUBROUTINE DIFF1 (RS, THS, RO, THO, ZO, AK, P , PI, D, N) 

2 COMPLEX PI, PR, PS, P, AI, AX 

3 AI - CMP LX ( 0 . , 1 . ) 

4 Cl « COS (.5* (THS - THO)) 

5 C2 = COS (.5* (THS + THO)) 

6 SI - SIN (.5* (THS - THO)) 

7 S2 = SIN (.5* (THS + THO)) 

8 RMR - (RS - RO) **2 + ZO**2 

9 R1K - SQRT (RMR + 4 . *Sl**2*RS*RO) *AK 

10 S1K - SQRT (RMR + 4 . *S2**2*RS*RO) *AK 

11 SQ1 - 2 . *SQRT (RS*RO) *AK 

12 IF (R1K .LT. l.E-5) THEN 

13 ARG1 - SQ1 *C1*1 . E5 

14 ELSE 

15 ARG1 - SQ1 *C1/R1K 

16 END IF 

17 IF (S1K .LT. l.E-5) THEN 

18 ARG2 - SQ1*C2*1.E5 

19 ELSE 

20 ARG2 - SQ1*C2/S1K 

21 END IF 

22 UR - ALOG (ABS (ARG1) + SQRT(ARG1**2 + 1 . ) ) *SIGN ( 1 . , ARG1) 

23 US - ALOG (ABS (ARG2) + SQRT(ARG2**2 + 1.) ) *SIGN ( 1 . , ARG2) 

24 CALL INTH (UR, R1K, D, N, PR, Jl) 

25 CALL INTH (US, S1K, D, N, PS, J2) 

26 AX - -AI*R1K 

27 PI = CEXP (AX) /R1K 

28 P - PR + PS 

29 RETURN 

30 END 


Comments: This subroutine calculates R. S, p^ and p, as defined in equation (3) from the input values of the source and 
observer coordinates. These parameters are then used in the subroutine INTH to calculate the two integrals in equation 
(3) multiplied by (-i/2). These two integrals are then added to give the final result 

4-7 Sine and cosine of the angles used in equation (3). 

8-10 Calculation of k-R and k-S where R and S are defined in equation (3). 

1 1 Calculation of factor in p^ and p, in equation (3). 

12-16 Calculation of the argument in the square brackets in equation (3) in the definition of p R . If the 

denominator in line 1 5 is very small (which will be true if the observer is very near the source) then a very 
large value is used for the argument to avoid division by zero. 

17-21 Same as 12-16 except for p, rather than p*. The argument in line 20 will be small if the observer is near 

the image source. 

22-23 UR = ^ and US = p,. 

24-25 Calculation of the integrals in equation (3) times factor (- i/2). 

26-27 Calculation of free-field source pressure times factor 4^/(0,^). See line 28 comment 

28 Source pressure with halfplane present times factor 471/(0^). 
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Appendix D: Subroutine for Transforming Sweep and Flow 

8/3/84 CORD by R. K. Amiet 

1 SUBROUTINE CORD (XS, YS, ZS, XO, YO, ZO, G, B, RSX, TSX, ROX, TOX, ZOX, XE, ZE) 

2 CA = COS (G/57 .2958) 

3 SA = SIN (G/57. 2958) 

4 SQ = SQRT (SA*SA + B*B*CA*CA) 

5 CX - B*CA/SQ 

6 SX - SA/SQ 

7 XOX = XO*CX - B*ZO*SX 

8 YOX = YO*B 

9 ZOX - XO*SX + B*ZO*CX 

10 XSX - XS*CX - B*ZS*SX 

11 YSX = YS*B 

12 ZSX - XS*SX + B*ZS*CX 

13 RSX - SQRT (XSX*XSX + YSX*YSX) 

14 ROX - SQRT (XOX* XOX + YOX* YOX) 

15 TSX = ATAN2 (YSX, XSX) 

16 TOX - ATAN2 (YOX, XOX) 

17 IF (TSX .LT. 0) TSX - TSX + 6.28318531 

18 IF (TOX .LT. 0) TOX - TOX + 6.28318531 

19 ZET = ZOX - (ZOX - ZSX) *ROX/ (ROX + RSX) 

20 XE - ZET*SX 

21 ZE - ZET*CX/B 

22 ZOX - ZOX - ZSX 

23 RETURN 

24 END 


Comments: This subroutine takes the coordinates for a non zero flow situation with airfoil sweep and transforms to a 
coordinate system with zero flow and zero sweep. It also calculates the coordinates of the edge point at which the ray 
strikes, in real coordinates. 

2-3 Sine and cosine of the sweep angle G or y . 

4-6 Sine and cosine of the transformed angle y ’ defined by equation (25). 

7-12 Coordinates x,y,z transformed first to X,Y ,Z system as in equations (17), then rotated through angle y ' as 

in equations (27). 

13-16 Transformed coordinates put in terms of polar coordinates r and 0. 

17-18 Angles to source and observer must be 0 < 0 < 2n. Negative angles not permitted because of the barrier. 

19 The value of £, the edge crossing point in equation (40), is calculated in a coordinate system with zero flow 

and the z axis along the airfoil leading edge. 

20-21 The coordinates of £ in the real system are calculated. 

22 The program INTH for no flow assumes that the z coordinate of the source is zero. 
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Appendix E: Subroutine for Simpson's Rule Integration 

8/12/84 INTGC by R. K. Amiet 


1 


SUBROUTINE INTGC (A 

2 


COMPLEX D, SUM 

3 


AN - N 

4 


DEL = . 5 * (B-A) /AN 

5 


N2 - 2*N - 1 

6 


X - A 

7 


M - 2 

8 


MB - 2 

9 


CALL FC (A, C, SUM) 

10 


CALL FC (B, C, D) 

11 


SUM - SUM + D 

12 


DO 50 I— 1, N2 

13 


X - X + DEL 

14 


M ■ -M 

15 


MB - MB - M 

16 


BM - MB 

17 


CALL FC (X, C, D) 

18 


SUM * SUM + BM*D 

19 

50 

CONTINUE 

20 


SUM - SUM*DEL/3 . 

21 


RETURN 

22 


END 


Comments: This subroutine performs a simple Simpson's rule integration. 

1 Inputs are 

A = lower integral limit 
B = upper integral limit 
C = constant used in integrand 
N = # of integration points 
4 DEL = 1/2 of step size. 

9-10 Find the value of the integrand FC at the two limits A and B. 

12-19 Evaluation of integrand at intermediate points. 

20 Multiplication by step size divided by 6. 
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Appendix F: Subroutine for Evaluation of Fresnel Integrals 

8/12/84 FRES by R. K. Amiet 


1 

2 

3 


SUBROUTINE FRES (X, E) 
DIMENSION A (12) , B(12), 
COMPLEX AI, E 
AI - CMP LX (0 . , 1 . ) 

10 


C ( 12) , D (12) 


7 

A(l) - 

1.59576914 

8 

A(2) - 

-1.702E-6 

9 

A(3) - 

-6.808568854 

10 

A (4 ) - 

-.000576361 

11 

A (5) = 

6.920691902 

12 

A (6) - 

-.016898657 

13 

A(7) - 

-3.05048566 

14 

A (8) - 

-.075752419 

15 

A(9) = 

.850663781 

16 

A ( 10) 

« -.025639041 

17 

A(ll) ' 

- -.15023096 

18 

A (12) 

- .034404779 

19 

B (1) * 

-3.3E-8 

20 

B (2) *= 

4.255387524 

21 

B (3) - 

-9.281E-5 

22 

B (4) - 

-7.7800204 

23 

B (5) » 

-.009520895 

24 

B (6) = 

5.075161298 

25 

B (7) = 

-.138341947 

26 

B (8) - 

-1.363729124 

27 

B<9) - 

-.403349276 

28 

B (10) 

= .702222016 

29 

B (11) 

- -.216195929 

30 

B (12) 

- .019547031 

31 

CALL SRIS (A, Y, 12, SA] 

32 

CALL SRIS (B, Y, 12, SB] 

33 

E - EXP (AI*X) *SQRT (' 

34 

GO TO 

20 

35 

10 Y - 4. 

/X 

36 

C<1) - 

0. 

37 

C<2) - 

-.024933975 

38 

C (3) - 

3.936E-6 

39 

C(4) - 

.005770956 

40 

C (5) = 

.000689892 

41 

C(6) = 

-.009497136 

42 

C (7) - 

.011948809 

43 

C<8) = 

-.006748873 

44 

C (9) = 

.00024642 

45 

C (10) . 

= .002102967 

46 

C(ll) ■ 

- -.00121793 

47 

C (12) . 

- .000233939 

48 

D (1) - 

.19947114 

49 

D(2) - 

2.3E-8 

50 

D(3) - 

-.009351341 

51 

D<4) - 

2 . 3006E-5 

52 

D<5) = 

.004851466 

53 

D(6) = 

.001903218 

54 

D(7) - 

-.017122914 
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55 D ( 8) - .029064067 

56 D (9) « -.027928955 

57 D (10) - .016497308 

58 D (11) - -.005598515 

59 D (12) - .000838386 

60 CALL SRIS (C, Y, 12, SC) 

61 CALL SRIS (D, Y, 12, SD) 

62 E = . 5* (1 + AI) + EXP (AI*X) *SQRT (Y) * (SC - AI*SD) 

63 20 RETURN 

64 END 


Comments: This subroutine is for calculation of the Fresnel integrals. It is based on reference 12. 

5 The calculation is divided into small and large x approximations. The small x (x < 4) approximation is 

given in lines 6-34. The large x (x > 4) approximation is given in lines 35-62. 

7-30 Input constants for the small x polynomial approximation. 

31-32 Evaluations of die polynomials for the small x approximation. 

33 Final evaluation of die Fresnel integrals C 2 (x) + iS 2 (x) for an input x. 

36-59 Input constants for the large x polynomial approximation. 

60-61 Evaluations of the polynomials for the large x approximation. 

62 Final evaluation of the Fresnel integrals C 2 (x) and S 2 (x) for an input x. 


Appendix G: Subroutine for Calculation of Integrand in Equation (8) 

S/12/84 FC by R. K. Amict 

1 SUBROUTINE FC (X, C,F) 

2 COMPLEX F, AI 

3 AI - CMP LX ( 0 . , 1 . ) 

4 Y - C*COSH (X) 

5 CALL BSL(Y,AJ0, AJ1, Y0,Y1) 

6 F - AJ1 - AI*Y1 

7 RETURN 

8 END 


Comments: This subroutine calculates the value of the integrand in equation (8). 
The inputs are X corresponding to p. in equation (8) and C corresponding to kR. 
4 Argument of the Hankel function. 

6 F is the Hankel function of the second kind. 


24 



Appendix H: Subroutine for Calculation of the Bessel Functions 

1976 BESL by R. K, Amict 

1 SUBROUTINE BSL (X, JO , Jl , YO , Yl) 


2 


REAL 

J0, Jl 


3 


PI 

* 

3.1415927 


4 


Y - 

> X/3 . 


5 


IF 

(Y-l.) 10, 10, 20 


6 

10 

Z2 

ac 

Y**2 


7 


Z4 

« 

Z2**2 


8 


Z6 

- 

Z4*Z2 


9 


Z8 

- 

Z4**2 


10 


Z10 

i - Z8*Z2 


11 


Z12 

: = Z6**2 


12 


JO 


1.- 2 . 2499997*Z2 + 1.2656208*Z4 - , 

.3163866*Z6 

13 

1 



+ . 0444479*Z8 - ,0039444*Z10 + . 

00021*Z12 

14 


Y0 

= 

2 . * J0*ALOG (X/2 . ) /PI + ,36746691 + . 

, 60559366*Z2 

15 

1 



- . 74350384 *Z4 + ,25300117*Z6 - 

,04261214*Z8 

16 

1 



+ . 00427916*Z10 - ,00024846*Z12 


17 


Jl 

- 

(.5 - ,5624 9985*Z2 + ,21093573*Z4 - 

- . 03954289*Z6 

18 

1 



+ . 004433 19*Z8 - . 00031761*Z10 + 

. 00001109*Z12) *X 

19 


Yl 

- 

2 , * Jl*ALOG (X/2 . ) /PI + (-.6366198 + 

. 2212091*Z2 

20 

1 



+ 2 . 1682709*Z4 - 1.3164827*Z6 + 

,3123951*Z8 

21 

1 



- . 0400976*Z10 + . 002787 3 *Z 12) /X 


22 


GO 

TO 30 


23 

20 

Z2 

- 

Y**2 


24 


Z3 

« 

Z2*Y 


25 


Z4 

» 

Z2**2 


26 


Z5 

= 

Z4*Y 


27 


Z6 

= 

Z3**2 


28 


F0 

* 

.79788456 -,00000077/Y - . 00552740/Z2 - . 00009512/Z3 

29 

1 



+ . 00137237/Z4 - .00072805/Z5 + 

. 00014476/Z6 

30 


TO 

« 

X-. 78539816 - . 04166397/Y- . 00003954/Z2 + . 00262573/Z3 

31 

1 



- .00054125/Z4 - .00029333/Z5 + 

. 00013558/Z6 

32 


S - 

* SQRT(X) 


33 


JO 

- 

F0*COS (TO) /S 


34 


Y0 

- 

F0*SIN (TO) /S 


35 


FI 

» 

.79788456 +.00000156/Y + . 01659667/Z2 +.00017105/Z3 

36 

1 



- . 0024 9511/Z4 + . 001 13653/Z5 - 

. 00020033/Z6 

37 


T1 

» 

X -2.3561945 +.12499612/Y +.0000565/Z2-.00637879/Z3 

38 

1 



+ . 00074348/Z4 + .00079824/Z5 - 

• 00029166/Z6 

39 


Jl 


Fl*COS (Tl) /S 


40 


Yl 

* 

F1*SIN(T1) /S 


41 

30 

RETURN 


42 


END 




Commons: This subroutine calculates the Bessel functions J, t Y 0 , Y r 

5 Two x ranges: x < 3 and x > 3. Small x range: lines 6-22, and the high x range: lines 23-40. 

12-15 Equations (9.4.1) and (9.4.2) of reference 1 1, pp. 369. 

16-20 Equations (9.4.4) and (9.4.5) of reference 1 1, pp. 370. 

28-34 Equation (9.4.3) of reference 11, pp. 369. Polynomial approximations. 

35-40 Equation (9.4.6) of reference 1 1 , pp. 370. 


25 


Appendix I: Program for Calling the Subroutine INTH 

7/21/84 CINTH by R. K. Amiet 


1 


COMPLEX P 

2 


CHARACTER* 10 FNAME 

3 


WRITE (9,*) ' OUTPUT FILE NAME 

4 


READ (9,*) FNAME 

5 


OPEN ( 1 , F ILE-FNAME , STATUS- ' NEW ’ ) 

6 

10 

WRITE (9, 100) 

7 

100 

FORMAT (/IX, ’ INPUT RK, D, N, J (CLOSE FOR J>0)'/) 

8 


READ (9,*) RK, D, N, J 

9 


IF (J .GT. 0) GO TO 20 

10 


WRITE (9, 300) RK, D, N 

11 


WRITE (1,300) RK, D, N 

12 

300 

FORMAT (/IX, 2F7. 2, 14) 

13 


DO 50 1-1,50 

14 


SQT - SQRT(RK) 

15 


TR - - . 1*I/SQT 

16 


G - TR/SQRT <2 . ) 

17 


UR - 2 . *ALOG (G + SQRT (G*G + 1.)) 

18 


CALL INTH (UR, RK, D, N, P, J) 

19 


APRK2 - RK*SQRT (2 . * (TR*TR+1) * (TR*TR+2) ) *CABS (P) 

20 


TRRK - TR*SQT 

21 


WRITE (9, 400) TR, TRRK, APRK2, J 

22 


WRITE (1,400) TR, TRRK, APRK2, J 

23 

400 

FORMAT (IX, 3E13.5, 14) 

24 

50 

CONTINUE 

25 


GO TO 10 

26 

20 

CLOSE (1) 

27 


END 


Comments: The program steps through TR corresponding to t R in equation (10b). The step size is adjusted in line 1 5 
according to the size of RK in order to collapse the approximate calculations into a single curve independent of the size of 
RK. This is done by making the step size in the argument of the Fresnel integral in equation (11) independent of RK. 
The results should be plotted against TRRK calculated in line 20; TRRK is proportional to the step number as desired 
while TR is not. In line 19 all the factors except CABS(P) are for normalizing the output to be independent of k R for 
the asymptotic solution. 

2-5 For inputting filename. 

6-9 Prompts for inputs of R k, D. N and J; inputting J > 0, closes the output file and ends the session. 

10-12 Writes the above inputs to the screen and to the output file for verification. 

13- 24 Loop to step through t r , the integral limit 

14- 15 t r is normalized so that the step size of the argument of F* is the same for all values of Rk. 

16-17 m is calculated from t R using equation (9). 

18 The integral I in equation (8) is calculated and multiplied by (- i/2). This is the integral in the equation 
following equation (3.2) of reference 5. 

19 CABS(P) is the absolute value of the integral in line 18. All the remaining factors are for cancelling factors 

other than F* in equation (1 1). The factor RK appears explicitely as Rk in equation (1 1). The two factors 
with TR-TR come from substituting for R,/R using equation (10). This can also be seen by backtracking a 
step to equation (8). For the approximation, t R is substituted for T in the radical in the denominator. The 
result must be multiplied by a factor of 1/2 since the program is calculating the text integral I, multiplied by 
a factor (- i/2). Finally, the factor 2 W comes from the limit of F*(x) 2‘® as x 0. 

20 TRRK is proportional to the step number, independent of the value of R k, whereas TR will not be. Since 
the argument of F* in equation (1 1) is proportional only to the step number, the result is plotted against 
TRRK so that the plot will be independent of Rk. 

21-23 The results are output to the screen and to the output file. 

25 Go back to prompt for new inputs. 
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Appendix J: Subroutine for Calculating a Power Series 

8/12/84 SR1S by R. K. Amici 


1 

2 

3 

4 

5 

6 

7 50 

8 
9 


SUBROUTINE SRIS (A, X, N, SUM) 
DIMENSION A ( 15) 

SUM = A (N) 

DO 50 I«1,N-1 
J - N - I 
SUM - SUM*X + A ( J) 
CONTINUE 
RETURN 
END 


Comments: This subroutine calculates a power series in the variable x given the array of constants A, the variable x and 
the number of terms N. The program uses Homer's rule; see e.g., reference 13. 

3 This term will be multiplied by x a total of N - 1 times. 

4-7 The A(J) term will be multiplied by x a total number of J - 1 times. 
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Half-plane 
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Figure 1: Geometry of the diffraction problem. The half plane lies in the y = 0 plane. Point A is the 
projection of the observer location on the x,y plane. The source lies in the x.y plane The 
zero sweep case is shown. 




29 




Figure 3: Definition of the sweep angle y. 



90 



Figure 5a: Normalized directivity of a point source diffracted by a leading edge in a flow; observer 
directly above the retarded position of the edge; i.e., 0. = 90. 
r. = 1, r„ = 1000, k, = 10, 0. = 141.34, M = 0.8, y= z. = 0. 
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Figure 6b: Normalized directivity of a point source diffracted by a trailing edge in a flow; source 
directly above the present position of the edge. 

r. = 1, r, = 1000, k. = 10, 0. = 90, M = -0.8, y = z. = 0. 
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Figure 7a: Mach 
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Figure 7b: Mach 
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